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Abstract. Matter-wave interferometers utilizing different isotopes or chemical 
elements intrinsically have different sensitivities, and the analysis tools available 
until now are insufficient for accurately estimating the atomic phase difference 
under many experimental conditions. In this work, we describe and demonstrate 
two new methods for extracting the differential phase between dual-species atom 
interferometers for precise tests of the weak equivalence principle. The first method 
is a generalized Bayesian analysis, which uses knowledge of the system noise to 
estimate the differential phase based on a statistical model. The second method 
utilizes a mechanical accelerometer to reconstruct single-sensor interference fringes 
based on measurements of the vibration-induced phase. An improved ellipse-fitting 
algorithm is also implemented as a third method for comparison. These analysis 
tools are investigated using both numerical simulations and experimental data from 
simultaneous 87 Rb and 39 K interferometers, and both new techniques are shown to 
produce bias-free estimates of the differential phase. We also report observations 
of phase correlations between atom interferometers composed of different chemical 
species. This correlation enables us to reject common-mode vibration noise by a 
factor of 730, and to make preliminary tests of the weak equivalence principle with 
a sensitivity of 1.6 x 10 -6 per measurement with an interrogation time of T = 10 
ms. We study the level of vibration rejection by varying the temporal overlap between 
interferometers in a symmetric timing sequence. Finally, we discuss the limitations of 
the new analysis methods for future applications of differential atom interferometry. 
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1. Introduction 

Einstein’s equivalence principle (EEP) is a fundamental concept in physics that describes 
the exact correspondence between the gravitational and inertial mass of any object. It 
is a central assumption of the theory of General Relativity—which interprets gravity 
as a geometrical feature of space-time, and predicts identical accelerations for different 
objects in the same gravitational field. Precise tests of the EEP are of great interest in 
various fields of physics. For instance, some theories that attempt to unify gravity with 
the other fundamental forces predict a violation of this principle [1, 2]. The detection of 
such a violation could aid our understanding of dark energy in cosmology, and advance 
the search for physics beyond the Standard Model. In contrast, null results are also 
pivotal for putting bounds on model parameters contained in various extensions to 
General Relativity [3, 4, 5]. The equivalence principle is generally divided into three 
sub-principles that each must be satisfied for the EEP to hold [ 6 , 7]: the local Lorentz 
invariance, the local position invariance and the weak equivalence principle (WEP). In 
this article, we will focus on the latter. 

The WEP—otherwise known as the universality of free fall—states that a charge- 
free body will undergo an acceleration in a gravitational field that is independent of its 
internal structure or composition. Tests of the WEP generally involve measuring the 
relative acceleration between two different test bodies that are in free fall with the same 
gravitational field. The WEP is characterized by the Eotvos parameter, r/, given by 

op — a 2 A a 

^ = 2 — -= —, ( 1 ) 

CL\ ~\~ CL 2 CL 

where a\ and 02 are the accelerations of the two bodies, A a = a\ — <22 is the relative 
acceleration, and a = (ai + a 2 )/2 is the average acceleration. The WEP is satisfied if 
and only if A a = 0—implying that p = 0. 

The most precise tests of the WEP have been carried out with lunar laser ranging 
techniques [ 8 ], or using a rotating torsion balance [9, 10], which have both measured 77 
at the level of a few parts in 10 13 . Various Space missions to test the WEP at improved 
levels (10 ~ 15 or better) using other classical devices are presently in progress [11, 12, 13]. 
On a separate frontier, a number of groups have carried out tests between cold atoms 
[14, 15, 16, 17, 18] in an effort to probe the WEP at the quantum level. The majority 
of these tests have been conducted using matter-wave interferometers which, over the 
past few decades, have been extensively studied both theoretically and experimentally 
[19, 20, 21 , 22 ], Atom interferometers have been utilized as ultra-precise inertial sensors 
to measure, for example, the gravitational acceleration g [23, 24, 25, 26], the gravitational 
constant G [27, 28, 29], gravity gradients [25, 30, 31, 32, 33], gravitational field curvature 
[34], and rotations [35, 36, 37, 38]. A WEP test based on atom interferometry involves 
measuring the differential phase shift resulting from a relative acceleration between two 
species with different masses that are in free fall within the same gravitational field. 
This measurement is based on the same principle as gravity gradiometry, where the 
quantity of interest is the differential phase between test atoms of the same type but in 
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different spatial locations. The gradient of the gravitational field can be extracted from 
the differential phase between two sources, while higher derivatives of the field can be 
accessed if more than two sources are used. This technique was recently demonstrated 
to measure the curvature of the gravitational field, and has been proposed to detect 
gravitational waves and to study geophysical effects [39, 40, 41]. Presently, the state-of- 
the-art for WEP tests using matter-wave interferometry corresponds to an uncertainty 
of 3 x 10 8 [18]. A comparison between the gravitational acceleration measured by 
atoms and a macroscopic object ( i.e. a falling corner-cube) have also been carried out, 
and yield agreement at the level of 8r) ~ 6.5 x 10~ 9 [42], A handful of ground-based 
[43, 44, 45, 46] and micro-gravity-based [41, 47, 48, 49, 50] cold-atom experiments are 
currently underway that aim to greatly improve this precision. In addition, there have 
been a number of proposals for Space-based quantum tests of the WEP [7, 51, 52, 53, 54] 
that target accuracies at the level of 10 -15 . 

So far, most tests with cold atoms have used two isotopes of the same atomic 
element, e.g. 85 Rb and 87 Rb [14, 15, 18, 55], or 87 Sr and 88 Sr [16]. Although this 
class of test bodies has demonstrated a good level of common-mode noise rejection 
when performing differential phase measurements [15], it is intrinsically less sensitive 
to possible violations of the equivalence principle because the two atoms are relatively 
similar in mass and composition. Thus, it is interesting to perform these tests with 
two entirely different atomic elements. In this article, we will focus on the case of 8 'Rb 
and 39 K. These atoms exhibit a large difference in their number of nuclei—facilitating a 
mass ratio of Mr^/Mk ~ 2.2. Additionally, they have identical hyperfine spin structure, 
and similar excitation wavelengths (around 780 nm and 767 nm, respectively), which 
enables the use of the same laser technology and optics for cooling and interferometry. 
Dual-species interferometers of this type have the added advantage of being highly 
independent—that is, atomic sample properties such as the size and temperature, or 
interferometer parameters such as the interrogation time, Raman phase, and detuning, 
can be controlled independently. In contrast to dual-isotope setups where many of these 
parameters are coupled, this feature is ideal for studying a variety of systematic effects 
that will be important for future precision measurements [55]. For a more complete 
comparison of alkali atoms as candidates for WEP tests, see for example ref. [56]. 

One complication that arises with non-common elemental species is a difference 
in the scale factors, Sj ~ hjT?, between the interferometers. When the interrogation 
times Tj are the same, this difference originates from the effective wave vectors kj S 
of the interferometer beams used for atoms j = 1 and 2. Assuming that the WEP 
is true, the phase shift of the two interferometers due to a common acceleration a is 
Tj = Sja. Thus, a difference in the scale factors produces a relative phase shift between 
interferometers for the same acceleration: S(p s J s — (Si — S^a. For the case of 85 Rb 
and 8 'Rb, the scale factors can be made the same by a suitable choice of Raman laser 
detuning that guarantees /cf ff = [52], However, this is not generally possible for 

different chemical elements, and the systematic phase shift must be addressed in other 
ways. 
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Another issue related to having different scale factors regards the rejection of 
common-mode vibration noise between interferometers. From an analysis of the 
interferometer transfer functions (see Appendix C or refs. [52, 57, 58], for instance), one 
can show that perfect common-mode rejection requires four conditions to be satisfied: 
(i) the interferometers occur simultaneously with Tj = T 2 , such that they experience 
the same vibration noise, (ii) they have identical wave vectors, kf 1 = k | ff , they exhibit 
(iii) identical effective Rabi frequencies, = f^®, and (iv) identical pulse durations, 
T\ — t 2 . These conditions imply that if S\ ^ S 2 , the interferometers do not respond to 
common-mode noise with the same phase shift. 

The scale factors can be made the same by adjusting the interrogation times of the 
interferometers such that T\ = rT 2 , where r = y/fcf®//cf® [48]. This technique eliminates 
the systematic phase shift 5<f> s J s resulting from a constant acceleration, and improves the 
rejection of common vibration noise at frequencies < l/Tf, but it degrades the rejection 
efficiency at frequencies above ~ 1/Tf (see Appendix C). However, if the ratio r is very 
close to unity, as it is for some choices of atoms (r ~ 1.009 for 39 K and 8 'Rb), this option 
represents a good compromise between efficient noise rejection and reducing systematic 
effects. 

In this article, we describe and demonstrate three analysis methods for atom- 
interferometric WEP tests—including two new techniques that eliminate both 
aforementioned problems of systematic phase shifts and diminished common-mode 
rejection between coupled interferometers of different atomic species. The first of these 
two new methods is a generalized Bayesian analysis of the Lissajous curves formed 
by plotting the coupled sensor measurements parametrically. The second technique 
involves restoring the interferometer fringes by correlating with an auxiliary mechanical 
accelerometer. In this case, the phase shift for each species can be measured directly 
from the reconstructed fringes regardless of their scale factors or the degree of temporal 
overlap between the interferometers. Both of these new methods intrinsically account for 
different scale factors, and return unbiased estimates of the differential phase. Finally, 
to give a complete picture, we compare these techniques with an improved ellipse-fitting 
method recently developed by Szpak et al [59]. This numerical procedure yields an 
estimate of the differential phase shift with reduced bias compared to more commonly 
implemented algorithms in the presence of significant amounts of uncorrelated noise 
between sensors. 

In this work, we also report correlated phase measurements between simultaneous 
interferometers of different chemical species ( 39 K and 8 'Rb). When operated in an 
environment with significant levels of vibration noise, we demonstrate a common-mode 
vibration rejection factor of 7 ~ 730. These results represent a major step toward precise 
tests of the WEP with elements exhibiting vastly different masses. We also investigate 
the accuracy of the three aforementioned methods on experimental data obtained from 
the K-Rb interferometer. 

The article is organized as follows. Section 2 reviews some theoretical background 
concerning a WEP test with a dual-species interferometer. In sec. 3, we briefly describe 
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the three methods of extracting the differential phase. We give a brief description of the 
experimental setup for the K-Rb interferometer in sec. 4. We present our experimental 
results in sec. 5, and we give a discussion of the advantages and limitations of the 
new methods in sec. 6. Finally, we conclude in sec. 7. A detailed description of the 
three analysis methods, including extensive numerical tests of the generalized Bayesian 
estimator, can be found in the Appendices. 

2. Testing the WEP with two atomic species 

An atom-interferometric test of the WEP involves measuring the relative acceleration 
between two atoms of different mass. This can be done in one of two ways: (i) 
the absolute acceleration of each atom, cq and < 22 , can be individually measured and 
subtracted, or (ii) A a can be measured directly from the differential phase, (fid- In 
the ideal case, acceleration measurements are performed simultaneously in order to 
take advantage of correlated noise between sensors—reducing the total uncertainty in 
Aa. Since method (ii) involves a direct measurement of (fid, it intrinsically requires 
both simultaneity and phase correlation between atomic sensors to reject common-mode 
noise. Henceforth, two or more atom interferometers that satisfy these conditions are 
referred to as “coupled sensors”. Method (i) can be carried out regardless of these two 
constraints. In this section, we outline some theoretical background related to a WEP 
test with method (ii). 

Generally, the output from two coupled atomic sensors is described by the following 
sinusoids 


yi(a) = A 1 cos(Sia + (fii) + Bi, (2a) 

V 2 ( a ) = A 2 cos(S 2 a + <f> 2 ) + B 2 , (2b) 

where Aj and Bj are, respectively, the amplitude and offset of the interferometer fringes 
associated with sensor j (j = 1, 2). In principle, these two parameters can be measured 
and eqs. (2) can be recast in the normalized form nj = {jjj — Bf)/Af 

'Ui(a) = cos(*Sia + (fi \), (3a) 

n 2 (a) = cos(S 2 a + (fi 2 ). (3b) 

Here, a is an acceleration common to both atoms, Sj is the scale factor for interferometer 
j, and (fij is a phase shift. The scale factors can be computed exactly from the integral 
of the response function, fj(t), given by eq. (C.5): 

Si = kf J fj(t)dt = if (T) + 275 ) (r, + tCj , (4) 

where kj S is the effective wave-vector for the counter-propagating interferometer beams, 
Tj is the interrogation time, and Tj is the 7 t /2 Raman pulse duration. A detailed 
explanation of the response function and its role in WEP tests is outlined in Appendix 
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C. For large interrogation times, Tj Tj, the scale factors reduce to the well-known 
relation Sj ~ k^Tj. 

The phases <f>j each have three contributions, one from the interferometer laser 
phase 0^ ser , one from parasitic systematic effects 0® ys , and one from the WEP signal 

0f EP 

*i=<i>'r+<i>T+*T SP - ( 5 ) 

In an experiment, the laser phase is a control parameter which can be set to zero, and 
systematic effects are independently nullified as much as possible. The shift due to a 
WEP violation can be defined as 0) VEP = Sj(aj — a), which is expected to be very close 
to zero. In the ideal case, the total interferometer phase contains only the shift due 
to the mean acceleration, Sja , and a WEP violation. It then follows that 

= Sjti + (f)j = SjOj. (6) 

In this case, the total phase of each interferometer can be related to the Eotvos parameter 
in the following way 

$i/S 1 -$ 2 /S 2 A a 

V = -= —• ( 7 ) 

a a 

In principle, the sensitivity in this type of WEP test increases as the square of the 
interrogation time, T ~ Tf ~ T 2 , due to the scale factors, Sj, that appear inversely in 
eq. (7). 

The general form of eqs. (3) describes a Lissajous curve. For the purposes of this 
analysis, it is useful to redefine the phases in eqs. (3) to reduce the number of free 
parameters. Choosing sensor 2 as a reference to rescale the phase of sensor 1, we define 
a common phase 0 C that satisfies 


4> c = S 2 a + 02 ) 
K,(j) c + 0d = Sid + 01- 


(8a) 

(8b) 


Here, we introduce two new parameters—the scale factor ratio k and the differential 
phase 0^—which are constrained from eq. (8) to be 


= 01 - K0 2 , 


K = 


Si 

S.f 


( 9 ) 


The sensor outputs are now recast with 0 C as the primarily parameter 


Til (0c) = cos (fi)0 c + 0d), (10a) 

n 2 (0 c ) = cos(0c). (10b) 


Comparing eqs. (6), (7) and (9), it follows that the Eotvos parameter is directly 
proportional to the differential phase: 


V = 


4>d 
~S\CL * 


( 11 ) 
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3. Correlative methods of differential phase extraction 

In this section, we review three different methods to measure the differential phase 
from experimental data: ellipse fitting, Bayesian analysis and fringe reconstruction from 
mirror acceleration measurements. 

3.1. Improved ellipse fitting 

The ellipse fitting technique was first applied to atom interferometry in ref. [60] for 
situations in which the phase common to two coupled atomic sensors is sufficiently 
scrambled to impede individual fringe observation. In this case, when the measurements 
from each sensor are plotted parametrically, one obtains an ellipse that is free from 
common phase noise. Using a least-squares ellipse fitting algorithm, the differential 
phase 4>d can be extracted. Multiple groups have demonstrated the utility of ellipse 
fitting for measurements of gravity gradients [32, 33] and the gravitational constant G 
[27, 28, 61]. However, this technique suffers from a number of drawbacks. First, it is 
valid only for coupled sensors with the same scale factor (k = 1). Second, in the presence 
of moderate amounts of noise in the fringe offsets or amplitudes [the parameters Aj and 
Bj in eqs. (2)], or in the differential phase, the ellipse fit returns a biased estimate of <f d 

t- 

Recently, Szpak et al [59] developed an algorithm based on the optimization of 
the approximate maximum likelihood distance which seeks a balance between costly 
geometric methods and stable algebraic techniques. This algorithm—termed the “fast 
guaranteed ellipse fitting” (FGEF) method—exhibits a smaller bias in the differential 
phase estimate over a relatively large phase range (centered on tt/2) compared to the 
more commonly used “direct ellipse fit” (DEF) technique [62]. Additionally, ref. [63] 
includes error estimations for the geometrically meaningful ellipse parameters (center 
coordinates, axes and orientation). We have extended their work to include an estimate 
of the statistical uncertainty in the differential phase, 5<f> d . We provide a more detailed 
comparison between DEF and FGEF methods of ellipse fitting in Appendix A. 

3.2. Generalized Bayesian analysis 

Heuristic approaches to estimating the differential phase, such as ellipse-fitting methods, 
do not have knowledge of the noise present in experimental data, nor of how various 
types of noise can affect the outcome of measurements. Bayesian analysis offers an 
efficient alternative to the problem by constraining the estimate based on a statistical 
model that describes the distribution of data that results from different noise sources 
[64], Bayesian phase estimation was studied in the context of atom interferometry in 
ref. [65] for two sensors containing the same scale factor (k = 1). In that work, a 

$ Rosi et al [34] demonstrated that the bias in the estimate of (f>d can be eliminated under 
certain conditions when fitting an ellipse in three dimensions from the output of three simultaneous 
interferometers. 
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Figure 1. The process of Bayesian estimation of the differential phase fd from 
synthetic data following a Lissajous curve, (a) Randomly chosen points (labelled 1 — 4) 
following eqs. (10) with Gaussian noise added to n\ and n 2 . For the actual Lissajous 
curve (shown as the black curve), we chose k = 0.8 and 4>^ ct = 1 rad for illustrative 
purposes, (b) The prior probability distribution computed from Bayes’ algorithm after 
each measurement. The vertical solid line indicates the differential phase used in the 
simulation, (c) Error in the Bayesian estimate after successive measurements. Points 
represent the difference between and ( i.e. the systematic error), where (jff 1 is 

the Bayesian estimate based on the maximum likelihood value from the corresponding 
prior distribution. The error bars indicate the statistical uncertainty, which are 
computed from the standard deviation of the prior distributions shown in (b). 


detailed study of each possible noise source (amplitude, offset and differential phase) is 
presented. Reference [56] also used Bayesian analysis to estimate the differential phase 
from a hypothetical system with k < 1. There, however, only noise in the differential 
phase is considered, and the range of common phase was constrained to 0 C G [0, 7 r]. To 
the best of our knowledge, no complete Bayesian estimator exists that (i) is valid for 
any scale factor ratio, (ii) accounts for noise in all relevant system parameters, and (iii) 
allows <f c to vary over a broad range. Furthermore, this type of analysis has not yet 
been demonstrated on experimental data from dual-species interferometers. 

In this work, we have developed a generalized Bayesian estimator for cfd —based on 
the approach of ref. [65]— that satisfies all three of the requirements mentioned above. 
We demonstrate this technique by measuring from both simulated data (see Appendix 
B) and experimental data from our K-Rb interferometer (see sec. 5). The advantage of 
using this estimation technique is that the uncertainty in cfd converges much faster than 
other methods (i.e. it scales as ~ l/\ f~N, where N is the number of measurements), so 
fewer data are required to reach a given level of sensitivity. Furthermore, since k is built 
directly into the Bayesian estimate of f>d , it is free from the aforementioned systematic 
phase shift 50[) ys arising between interferometers with different scale factors. However, 
some of the drawbacks of the Bayesian analysis are that it requires a priori knowledge 
of the noise in the system, and it is computationally costly due to the large number of 
integrals that must be evaluated. 

Figure 1 illustrates the basic Bayesian estimation procedure. Here, we simulate 
data that follow the Lissajous equations (10) with added Gaussian noise in the sensor 
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offsets. After each successive measurement from the system, the width of the new “prior” 
probability distribution decreases and additional peaks are suppressed—facilitating 
an improvement in the estimate of 0^. This is how the Bayesian method builds in 
information from previous measurements. It is clear from figure 1(c) that after only a few 
iterations, both the statistical and systematic error in 0^ have decreased dramatically. 
A detailed description of the generalized Bayesian analysis can be found in Appendix 


B. 


3.3. Fringe reconstruction by accelerometer correlation - The differential FRAC 
method 

Differential atom interferometry is often utilized under conditions where each sensor is 
overwhelmed by external phase noise that is common to both sensors. Typically, one is 
concerned with only the differential phase and not the common phase 0 C , which is treated 
as an arbitrary parameter. Both the ellipse-fitting and Bayesian estimation methods 
for extracting 0^ take this approach. An alternative technique involves measuring the 
common phase and correcting for it. For the case of parasitic mirror vibrations, single¬ 
sensor interference fringes that are otherwise smeared by phase noise can be restored 
based on measurements from seismometers [66, 67, 68] or mechanical accelerometers 
[41, 69, 70]. Henceforth, we refer to this as the fringe reconstruction by accelerometer 
correlation (FRAC) method. In this work, we demonstrate how the FRAC method 
can be applied to two quasi-simultaneous interferometers of different atomic species to 
measure the relative phase shift between them. This technique to extract 0^ is referred 
to as the differential FRAC method throughout the article to differentiate between the 
(standard) FRAC method, which is generally employed to measure the absolute phase 
shift of a single atom interferometer. 

Figure 2 illustrates the basic schematic of the FRAC method for a single 
interferometer. A mechanical accelerometer is secured to the back of the reference mirror 
used to retro-reflect interferometry light, and the time-dependent mirror acceleration, 
a vlh (t), is recorded during the interferometer sequence. These acceleration measurements 
are first weighted by the response function of the j th interferometer, fj(t), and are then 
integrated to find the vibration-induced phase given by§ 



( 12 ) 


For each repetition of the experiment, this random phase is computed and correlated 
with the interferometer signal. This process allows one to reconstruct the interference 


fringes point-by-point. Depending on the level of vibrations and the interferometer 


§ The underlying assumption of this technique is that motion of the reference mirror at frequencies 
within the response bandwidth of the mechanical accelerometer are responsible for phase changes of 
the atom interferometer. Although the corresponding acceleration signal is indistinguishable from 
fluctuations in local gravity (as a consequence of the equivalence principle), we can be confident in our 
assumption since typical variations in gravity occur on timescales much larger than 2 T. 
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Acceleration Measurements 



Figure 2. Schematic of the FRAC method for a single interferometer. A mechanical 
accelerometer is mounted to the back of the retro-reflection mirror for the Raman beam 
(with wave vectors k\ and /C 2 ). Acceleration measurements during the interferometer 
sequence are weighted by the response function /(f) [see eq. (C.5)] and integrated to 
obtain the phase estimate /> vlb . Correlating this phase with the interferometer signal 
during the same time interval reproduces the interference fringe. 


sensitivity, the range of vibration-induced phases can span multiple fringes—enabling 
the single-sensor phase shift (pj to be measured using, for instance, a sinusoidal least- 
squares fit to the data, ft is straightforward to extend this algorithm for two or more 
interferometers, which do not need to be overlapped in time. In this case, the only 
additional requirement is that the time-series of mirror acceleration measurements span 
the interrogation times for all interferometers. For two coupled sensors, the differential 
phase is easily computed from the individual sensor phase shifts via (pd = <f> 1 — The 
statistical error in this quantity is governed by 

(S4>d) 2 = {Sf> 1) 2 + (nScj) 2 ) 2 - 2np (j , 1 ^ 2 (6(j) 1 )(S(j)2), (13) 

where the 8(pj represent the statistical uncertainties in the (pj obtained from fits to the 
two fringes, and is the correlation coefficient for the measurements of (p\ and </> 2 . 

In the limit of perfect correlation = 1), the uncertainty in the differential phase 

reduces to 8<pd = \8<pi — nSfol- Figure 3(a) illustrates how the coupled-interferometer 
correlation is utilized by the differential FRAC method. Since the fringes for each 
interferometer are recovered using measurements from the same classical device, the 
phase noise present on each fringe is highly correlated. This induces a correlation 
between the measurements of (p± and 0 2 extracted from the fits, as characterized by 
£? 0 i, 02 - The key to the differential FRAC method is maximizing this correlation to 
reduce the uncertainty in <f>d■ 
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Figure 3. (a) Schematic illustrating the source of correlation in the differential FRAC 
method. The signal from each atom interferometer (AI), exhibiting a known scale 
factor Sj and unknown phase shift fij , is directly correlated (indicated by the solid 
lines) with a common vibration phase measured by a mechanical accelerometer (MA). 
Since the resulting fringes are derived from a common source, measurements of each fij 
are highly correlated, as indicated by the dotted lines and characterized by the positive 
coefficient This technique can be extended for multiple coupled interferometers 

(shown in gray), although we focus on the case of only two. (b) Correlation coefficient 
Q<t>i, 4>2 as a function of fid estimated from a large sample of simulated fringes. Here, 
we assumed k = 1 and we added non-common Gaussian noise to the phase of the 
reconstructed fringes with a standard deviation of a ( j >d =0.1 rad. 


The correlation coefficient for a given set of reconstructed fringes can be estimated 
numerically from a large sample of simulated data. We find that it is sensitive to 
experimental parameters such as the level of uncorrelated noise on each sensor, the scale 
factor ratio and the differential phase. For instance, figure 3(b) shows the dependance 
of o ( j n )( j >2 on fid for synthetic fringes that contain non-common phase noise with a 
standard deviation of 0.1 rad. The correlation coefficient yields a maximum when the 
interferometers are perfectly in-phase or 7r radians out-of-phase. This is an ideal feature 
for WEP tests, since the maximum sensitivity occurs exactly at the expected signal 
of fid = 0. This implies that, unlike ellipse-fitting methods where the sensitivity is 
optimized at fid = 7t/2, one does not need to engineer an additional phase shift between 
the atoms to optimize the sensitivity and reduce systematic bias. Furthermore, a recent 
study of a gradiometer configuration (i.e. k — 1) has shown that the differential FRAC 
method can reach sensitivities close to the quantum-projection-noise limit when modest 
levels of uncorrelated phase noise are present [71]. 

A number of ideal features make this technique interesting for both absolute and 
differential atom interferometry experiments. 

1) The differential FRAC estimate of fid is precise and unbiased over the full phase 
range fid £ [0,7r], since it relies on least-squares fits to individual fringes. 

2) It is simple, fast, and computationally low in cost—allowing the interferometer 
phase to be corrected in real-time [70], or by post-processing the data [66, 67, 69]. 

3) Unlike the Bayesian analysis, the FRAC method does not require any a priori 
information about the interferometer offsets, contrasts, and noise parameters— 















Correlative methods for quantum tests of the weak equivalence principle 


12 


which can be challenging to measure accurately without phase stability [72]. 

4) Systematic phase shifts in <f>d due to non-identical pulse durations Tj and Rabi 
frequencies [15] are accounted for in the estimates of <f>f h for each interferometer. 
Such systematics will be important to consider in future long-baseline differential 
interferometry experiments [7, 44, 45, 46, 51, 52, 53, 54], 

5) The relative timing between coupled interferometers can be freely chosen—they 
need not be overlapped. This is a unique feature to dual-species interferometers 
that do not share the same Raman beams. Unlike the ellipse-fitting and Bayesian 
techniques, the FRAC method allows one to extract absolute phase information 
from each sensor. Varying the temporal overlap between interferometers can be 
useful for studying a variety of effects, such as the level of correlation between 
sensors, or systematics related to the interaction between atoms [55]. 

6 ) Single-sensor fringes can be accurately measured in “noisy” environments, which is 
ideal for mobile sensors such as atomic gravimeters [26, 41, 68, 69]. 

Although the standard FRAC method is conceptually simple to implement, the 
drawback is that it is sensitive to errors in the measurements of vibrations. Such 
errors include the quality of coupling between the mirror and the mechanical device, 
electronic noise in the signal acquisition, the level of self-noise of the device, drifts in 
the offset or sensitivity factor, and non-linearities in both the amplitude and frequency 
response. The natural low-pass filtering feature of atom interferometers can alleviate 
some of these effects—particularly those that dominate at frequencies beyond the cut-off 
frequency 1/2 Tj. Alternatively, the mechanical accelerometer can be used as a course 
measurement of the phase to identify the correct interferometer fringe. Then one can use 
the atom interferometer output to refine the phase measurement [69, 70] by inverting the 
sinusoidal relation (2). On the other hand, measurements of <fd using the differential 
FRAC method are much less sensitive to many of these noise sources since they are 
common to two simultaneous interferometers. We discuss the limitations of this method 
in more detail in sec. 6. 

4. Description of the ICE experiment 

ICE (Interferometrie Coherente pour l’Espace) is an experiment that aims to measure p 
using a dual-species interferometer of 8 'Rb and 39 K. It is designed to be transportable 
and to operate in the micro-gravity environment provided by the Novespace Zero -g plane 
[41, 48, 69, 73]. In this section, we give a brief description of the experimental setup. 

A detailed description of the telecom-frequency fiber-based laser system used on 
ICE can be found in refs. [41, 74], For each atomic species, we utilize a master-slave 
architecture, where the master laser diode is locked to either a saturated absorption 
peak (in the case of rubidium), or to a frequency comb (in the case of potassium). The 
slave lasers are frequency-locked to their corresponding master through an optical beat- 
note in the 1550 nm telecom band. After second harmonic generation to 780 nm for 
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8 'Rb and 767 nm for 39 K, the frequency of each slave laser can be precisely adjusted 
over ~ 1.3 GHz within ~ 2 ms of settling time. Approximately 1.5 W of total light is 
available in each slave beam before entering a free-space optical bench. This module is 
composed of a series of shutters and acousto-optic modulators (AOMs) that are used 
to split, pulse and frequency shift the light appropriately for cooling, state preparation, 
interferometry and detection. Finally, the 780 and 767 nm light is coupled into a series of 
single-mode, polarization-maintaining fibers and sent to the vacuum chamber. The two 
frequencies required for cooling and repumping, as well as driving Raman transitions in 
87 Rb, are generated via a broadband fiber-based electro-optic modulator operating near 
6.8 GHz. Similarly, an AOM operating in dual-pass configuration at ~ 230 MHz is used 
to generate these frequencies for 39 K. 

The sensor head is composed of a non-magnetic titanium vacuum chamber 
surrounded by a /j-metal shield. The chamber resides within three nested Helmholtz coils 
used to compensate residual magnetic fields and to generate a bias along the vertical axis. 
A custom 2-to-6 way fiber splitter is used to combine the 780 and 767 nm light intended 
for laser cooling without significant power loss via a polarizing cube and a dichroic wave 
plate. The splitter subsequently divides the light equally into six beams that are re¬ 
coupled into independent fibers used for the dual-species vapor-loaded magneto-optical 
trap (MOT). In a similar way, light for both detection and interferometry is overlapped 
in a free-space 2-to-l way fiber combiner for 780 and 767 nm. The ~ 2 cm diameter 
beams output from the combiner have the same linear polarization, and are aligned 
along the vertical direction through the atoms. A quarter-wave plate (fabricated for the 
intermediate wavelength of 773 nm and mounted in front of the retro-reflection mirror) 
rotates the polarization of the Raman beams by 90° such that the counter-propagating 
fields have linTlin polarization. 

A typical experimental sequence for the K-Rb interferometer is shown in figure 4 
and is carried out as follows. The MOT beams load approximately 2 x 10 8 (7x10') atoms 
in 0.5 s, which is followed by a 7 ms (5 ms) molasses cooling stage for the 8 'Rb ( 39 K) 
sample. In addition to cooling, the rubidium molasses stage also pumps the atoms into 
the |F = 2) ground state. This is followed by a microwave 7r-pulse that transfers atoms 
into \F — 1, mp — 0), and the remaining atoms are removed with a push beam resonant 
with the F = 2 to F' = 3 transition. During the potassium molasses, the frequency 
and intensity of the cooling and repump beams are modified in a similar manner to 
refs. [75, 76]. At the end of the molasses, the atoms are in a superposition of both 
hyperhne ground states, which is a critical part of the cooling mechanism for potassium 
[75]. We detune our 767 nm push beam to the red of the F = 2 to F' = 3 transition 
by r\j 17 MHz (2.9 T) to optically pump the atoms into the F — 1 level with a 3 /rs 
pulse. Following this depumping stage, the 39 K atoms are distributed roughly equally 
amongst the magnetic sub-levels of the lower hyperhne ground state. With this system, 
we achieve temperatures of ~ 3 pK for 8 'Rb and ~ 20 pK for 39 K, as confirmed by 
both time-of-hight imaging and velocity-sensitive Raman spectroscopy. After preparing 
the internal atomic states, we typically wait ~ 12 ms for the atoms to fall such that 
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Figure 4. Symmetric timing sequence for the K-Rb interferometer. The first Raman 
pulse of the 39 K interferometer is delayed by TRb + ATRb.K relative to that of 87 Rb 
such that the central 7r-pulse of both interferometers occurs at the same time. The 
preparation stages (shown in blue) include both the internal state preparation pulses 
(microwave + push beam for 8| Rb, depump for 39 K), and the external state preparation 
via time-of-flight. Raman pulses are shown in green, and detection pulses in orange. 


the Doppler resonance of both sets of counter-propagating Raman beams becomes non¬ 
degenerate. Additionally, we apply an external magnetic bias field between 1 — 2 Gauss 
to shift the |F = 1, m f = ±1) states of potassium away from the central itif = 0 state on 
which we perform interferometry. The frequency of the Raman beams for both species 
is detuned by —1.2 GHz (—200 T) relative to the F = 2 to F' = 3 transition. We then 
apply the interferometry pulses in a symmetric fashion, such that the central 7r-pulse for 
both interferometers occurs at the same time, as shown in figure 4. The delay between 
the 7t/2 pulses for either atom, AT R b,K, can be adjusted within the interrogation time 
of the rubidium interferometer, T Rb , in order to study correlations and effects related to 
the scale factor ratio, k. Finally, we measure the atomic state populations for each atom 
via fluorescence detection on an avalanche photodiode (50 MHz bandwidth) within 100 
gs of one another. 

5. Experimental results 

We now describe some experimental results obtained from the K-Rb interferometer. All 
of the data presented in this work were recorded in a laboratory environment, with the 
interferometer beams aligned along the vertical direction, and with no anti-vibration 
platform. To compensate for the Doppler shift due to gravity, the frequency difference 
between Raman beams for interferometers j — 1 = K and j — 2 = Rb is chirped at a 
rate of olj ~ kj S g to account for the gravity-induced Doppler shift of the falling atoms. 
This modifies the total phase shift of the interferometers from eqs. (2), = Sja + c , 

to the following 


dq = Sj(a - otj/k'f) + fa ~ {k!f g - aj)T?. 


( 14 ) 
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2n (a - ao) T 2 + ^ vlb (rad) 

Figure 5. Interferometer fringes from the s 'Rb interferometer operating at T = 25 
ms at a contrast of ~ 40% without vibration isolation. The interferometer phase is 
scanned by both uncontrolled mirror motion, and by varying the chirp rate a between 
Raman beams about a central value of ao = 25.1355 MHz/s. The open circles indicate 
raw measurements of the normalized atomic population in the |F = 2) state, while the 
closed circles are the same measurements after applying the FRAC phase correction, 
f vlh , to each point. A few of these corrections are shown as blue arrows. The solid 
curve is a least-squares fit to the corrected data, resulting in a signal-to-noise ratio 
of ~ 30 and a relative statistical uncertainty of 10 _ ' in the determination of g Rb— 
corresponding to almost an order of magnitude improvement compared to the raw 
data. 


The last expression represents the case when both interferometers experience the same 
acceleration, aj = a = g, and the scale factors can be approximated as Sj ~ 
Determining the location of the central fringe, for which aq = kj S g is fixed for all Tj , 
yields a measurement of g. Using the fact that the sensitivity of the interferometer 
scales as Tjy absolute measurements of g have been demonstrated at accuracies of a few 
10~ 9 [24, 68, 77, 78], 

As discussed in the introduction, we are interested in measuring the differential 
acceleration A a between 39 K and 8 'Rb. One way of achieving this is to measure the 
gravitationally-induced accelerations g k and g Rb from each interferometer independently 
by scanning the chirp rates, oq, in a low-noise environment. This is the approach 
recently employed for WEP tests with 39 K and 87 Rb by Schlippert et al [17]. However, 
at high levels of sensitivity (i.e. large Tj), or in “noisy” environments, mirror vibrations 
can corrupt the fringes—making individual phase measurements more challenging. We 
now demonstrate the utility of the FRAC technique for measuring g from a single 
interferometer under these conditions. 

There are typically two approaches in which the FRAC method can be applied to 
restore the interference fringes of a single interferometer. The first approach is to let 
the interferometer phase be “scanned” randomly by vibrations while the laser-induced 
phase is held fixed. The reconstructed fringes in this case are purely a function of 
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</>J lb , as shown in figure 2. This mode of operation can be used to precisely calibrate 
the mechanical accelerometer by rescaling the voltage-to-acceleration sensitivity factor 
of the device such that the fringe period is 27r||. The second approach is to scan 
the interferometer phase in a controlled manner, for example by varying the phase 
difference between Raman lasers, and to correct each phase using 0 vlb obtained during 
the same measurement interval. This procedure is illustrated in figure 5, where the 
fringes of a T = 25 ms 8 'Rb interferometer are shown before and after applying the 
FR.AC correction. Here, the interferometer is operated without any vibration isolation 
in the presence of a root-mean-squared (rms) DC vibration noise of a/* s ~ 6 x 1CT 5 
m/s 2 (integrated over the frequency response of the interferometer)—corresponding to 
an rms phase noise of = kf®/if^ s Tf ih ~ 0.6 rad. Acceleration measurements were 
performed with a force-balance three-axis accelerometer (Nanometrics Titan, DC to 
430 Hz bandwidth, 5 V/g sensitivity). By applying the FRAC correction to these data, 
we improve the signal-to-noise ratio (SNR) and hence the uncertainty in the central 
fringe measurement by almost an order of magnitude. We estimate an individual phase 
correction uncertainty of 5<p v ' h = 1/SNR ~ 33 mrad based on the improved SNR of 
~ 30. With this method, we emphasize that the interferometer sensitivity is directly 
linked to the intrinsic noise of the accelerometer + signal acquisition system, and the 
quality of the coupling between the device and the Raman mirror. Therefore, modest 
improvements to any of these system components can result in a dramatic increase in 
the fringe SNR. 

5.1. K-Rb Interferometer Correlation 

Typically, when mirror motion is the dominant source of phase noise it is advantageous 
to use differential atom interferometry techniques to measure A a through the differential 
phase <fd- This requires a high level of correlation between interferometers in order to 
reject the common-mode phase noise. We now compare three methods of extracting 
(fd from experimental data recorded in an environment with high vibrational noise, 
as in the case of onboard applications [58, 69]. These studies are also applicable to 
future high-sensitivity differential interferometers operated in low-noise environments 
[44, 45, 46], 

Figure 6 shows data produced by quasi-simultaneous K-Rb interferometers at a 
total interrogation time of 2 T = 6 ms. Here, we held the chirp rate fixed at otj ~ k^g 
for each species, and we applied strong vibrations to the system (a/* s ~ 0.05 m/s 2 ) 
such that the random vibration-induced phase 0J lb spanned multiple fringes (0/*s — 7.3 
rad). The vibrations were applied by mounting a heavy industrial fan on top of the 
support structure surrounding the vacuum system and running it during the experiment. 
Figure 6(a) shows a histogram of 8 'Rb |F = 2) population measurements, t/Rb, which 

|| One advantage of performing this procedure is that the device can be precisely calibrated for the 
vibration spectrum on site. Depending on the bandwidth and spectral response of the device, the 
sensitivity can vary significantly with the vibration spectrum. 
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Figure 6. Measurements of normalized \F = 2) state populations from simultaneous 
K-Rb interferometers operating at T ~ 3 ms. Fringes are scanned randomly by 
applying vibrational noise to the retro-reflection mirror. Graphs (a) and (b) show a 
time series of 500 measurements from the 8 'Rb and 39 K interferometers, respectively, 
along with histograms of the populations, (c) Atomic populations from (a) and (b) 
plotted parametrically—indicating strong correlation between the two species. The 
solid green line is an ellipse fit to the data using the FGEF method, which yields a 
differential phase cj>ff lpse = 1.13(2) rad. A separate estimate from a Bayesian analysis 
gives 0 Rayes = 1.18(2) rad. (d) Interferometer fringes reconstructed from measurements 
of mirror motion using the FRAC method. The red and blue curves correspond to least- 
squares fits to Rb and K data, respectively. The differential phase estimated from 
the fits is <(> RRAC = 1.17(1) rad. Other interferometer parameters: pulse separations: 
Tr b = 3.018 ms, Tk = 3 ms; 7r/2-pulse durations: TRb = 4 /is, tr = 6 ps; delay between 
interferometers: AXk,Rb = 10 ps; one-photon Raman detunings: ARb = Ar — —1.2 
GHz. 


clearly indicates the characteristic bimodal probability distribution of a sinusoid. These 
distributions can be used to estimate the contrast, offset and SNR of the interferometer 
fringes as described in ref. [69]. We note that the bimodal distribution is less pronounced 
for 39 K in figure 6(b) owing to a smaller fringe contrast, and thus a lower SNR, compared 
to s 'Rb. Despite this fact, the two sensors exhibit strong correlations, as confirmed by 
the ellipse in figure 6(c). 

For these experimental parameters the scale factor ratio is n = Sk/S Rb = 1-008, 
and the Lissajous curve formed by parametrically plotting the atomic state populations, 
?/Rb and Dk, is indistinguishable from an ellipse at the present level of offset noise. We 
measure a differential phase of 0^ lhpsc = 1.13(2) rad from a least-squares fit to an ellipse 
using the FGEF method [59]. We also estimate 0^ ayCh = 1.18(2) rad using the Bayesian 
analysis described in sec. 3.2 and Appendix B. Here, it is worth mentioning that this 
non-zero differential phase does not originate from a WEP violation, but from systematic 
phase shifts in the experiment—primarily due to the quadratic Zeeman effect from an 
external magnetic bias field (~ 1 G) that is used to sufficiently split the ground state 
magnetic sub-levels in 39 K. 

Figure 6(d) shows the output of each interferometer as a function of the vibration- 
induced phase, </>J lb . Here, the single-sensor fringes were reconstructed using the 
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FRAC method using mirror vibration measurements from a broadband micro-electro¬ 
mechanical accelerometer (Colibrys SF3600, DC to 1 kHz bandwidth, 1.2 V /g 
sensitivity). From these data the differential phase shift between interferometers is 
clearly visible. Sinusoidal least-squares fits to each fringe yield 0™ AC = 4>k — «0Rb = 
1.17(1) rad. Here, the statistical uncertainty 50 RRAC was computed from the quadrature 
sum of each interferometer phase error. The value of <f>d estimated from the Bayesian 
analysis and the FRAC method are in good agreement. On the other hand, the 
differential phase from the ellipse fit is underestimated by ~ 40 mrad, i.e. 2 a below 
^Baych 0 FRAC. yy e attribute this discrepancy to the inherent bias of ellipse-fitting 
techniques (see Appendix A), which increases with the level of offset noise or differential 
phase noise in either interferometer. 

We emphasize that a crucial input parameter for the Bayesian analysis is the 
common phase range. We use the accelerometer data to estimate this range once the 
experiment is complete: 0 C G [min(0^),max(0^)]. However, if an accelerometer is 
not available, it is also possible to estimate this range using the raw data from a single 
interferometer. For example, one can reduce the interrogation time until the sensitivity 
to vibrations reaches a point where interference fringes are clearly visible. By measuring 
the rms scatter of the phase about a reference sinusoid, one can estimate the level of 
vibration noise via the relation a/* s = (p^ s /S. Once is known, this relation can 
be inverted to determine the range of phase scanned by the same level of vibrations at 
larger sensitivities/interrogation times. 

The data shown in figure 6 (d) also indicate that the combined differential-atomic¬ 
sensor + mechanical-accelerometer system is capable of efficiently rejecting common 
vibrational noise. We estimate a rejection factor of 7 = /c efr a))*T 2 /(50^ RAC ~ 730 for 
these data. 

Figure 7 displays the results of a correlation study between rubidium and potassium 
interferometers operating at a total interrogation time of 2 T = 20 ms. Similar to figure 6 , 
the interferometer phases are scanned by externally applied vibrations (a)'* s ~ 1.7 x 10 -3 
m/s 2 , ~ 2.7 rad at T Rb = T K = 10 ms). The vibrations were applied using the 

same method as previously mentioned, but with the fan set on a slower rotation setting. 
Here, we vary the interrogation time of potassium, T K . in a symmetric way with respect 
to rubidium such that the centers of the 7 r-pulses coincide. This optimizes the degree to 
which the vibration-induced phase noise remains common-mode, while modifying the 
degree of temporal overlap between interferometers. It also allows us to control the scale 
factor ratio since k scales at (Tr/Tri,) 2 - 

From figure 7, three features are clearly visible as Tr is decreased. First, the 
potassium fringes undergo a phase shift that modifies the differential phase relative to 
the rubidium fringes. This feature, along with the fact that the scale factor ratio is 
varied, causes the shape of the Lissajous figures to change, as shown by the solid green 
curves. Second, the phase range scanned by the potassium interferometer reduces, since 
it scales as Tr. Finally, the level of correlation between the interferometers degrades as 
the temporal overlap decreases. This is evident from the lack of agreement between the 


Correlative methods for quantum tests of the weak equivalence principle 


19 


(a) 


T Rb = 10.09 ms, T k = 10 ms, ?c = 1.00 




(C) 


« « b (rad) 

T Rb = 10 ms, T K = 9 ms, k = 0.82 



-6 -4-2 0 2 

j-vib j,yib . 



(e) 


(rad) }W> 

T Rb =10ms, T K = 7ms, k = 0.50 



-6 -4-2 0 2 

j-vib j.vib . 


« «' 6 (rad) 



(b) 


T Rb =10ms, T c = 9.5 ms, k — 0.92 



(d) 



« ib (rad) 

T Rb =10ms, T K = 8ms, k = 0.65 




(f) 


« b (rad) VRb 

T Rb =10ms, T K =6ms, k=0.37 



-6-4-20 2 4 6 

< b «! b (rad) 



Figure 7. Correlation studies for different levels of temporal overlap. The symmetric, 
quasi-simultaneous K-Rb interferometer was operated with Tr b = 10 ms and the 
interrogation time for potassium was varied between Tk = 6 — 10 ms. The 
interferometer phase was scanned by externally applied vibrations and individual 
fringes were restored using the FRAC method. Parametric plots of the atomic 
populations are shown to the right, along with the expected Lissajous curve (solid green 
line). These curves result from plotting the fit functions to each reconstructed fringe 
parametrically. There is a clear disagreement between the predicted Lissajous curves 
and the data for Tr b — Tk > 2 ms. Other interferometer parameters: TRb = tk = 3 
p,s; ARb = Ar — —1.2 GHz. 


data and the predicted Lissajous curves, particularly for T K < 8 ms. 

Regardless of this degradation of correlation and temporal overlap between 
interferometers, the differential FRAC method is able to restore the interference fringes 
with a good SNR (~ 30 for 8 'Rb, ~ 10 for 39 K, limited by uncorrelated offset noise). 
The vibration rejection factor for each of the data sets shown in figure 7 is approximately 
7 ~ 100. This permits unbiased estimates of 0^ with a statistical uncertainty at the 
level of 5<f>d ~ 25 mrad with 300 points—corresponding to WEP test with a statistical 
sensitivity of 5r] ~ Scfd/k^gT^ = 1.6 x 10 -6 per data set. The robustness of the 
differential FRAC technique under “noisy” conditions makes it an ideal candidate for 
future WEP tests [7, 53, 54], or other differential atom interferometry applications 
[34, 40], 
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Figure 8. Comparison of differential phase estimates from the Bayesian and FRAC 
methods for various fd and u. (a) Difference between Bayesian and FRAC estimates 
^Bayes_^FRAC ag a f unc ti 0 n 0f RAC = 0 k —K</>Rb- The vertical error bars represent the 
combined statistical uncertainty of both estimates, (b) Difference between Bayesian 
and FRAC estimates as a function of k = Sk/SjO) corresponding to each point in (a). 
The majority of the data points shown in both figures are consistent with zero to within 
lcr of uncertainty, indicating agreement between the two techniques at the level of ~ 40 
mrad. The standard deviation of offset noise for each data set was typically ctb k ~ 0.20 
and cr s Rb — 0.05 in the normalized space ( n\ and 122 ). A value of a^, d = 0.05 rad was 
used for the differential phase noise of all data sets. The range of common phase noise 
was estimated from accelerometer measurements. 


In contrast, for Bayesian estimation, an increase in uncorrelated phase noise is 
problematic. When the “common phase” becomes largely uncorrelated, the Bayesian 
method can converge on multiple possible f>d, or may not converge at all. For these data, 
we find that by T K = 9 ms the Bayesian estimate of (fid is not consistent with the FRAC 
estimate, and for Tk < 8 ms the analysis is not able to converge on a unique value. We 
note that these particular results are strongly dependent on the level of phase noise, the 
degree of temporal overlap, the value of <fid and the scale factor of each interferometer. 
In the following section, we study some of these dependencies more quantitatively. 

5.2. Comparison of Bayesian and FRAC methods as a function of n and (fid 

We have tested the functionality and accuracy of both the Bayesian and FRAC methods 
for extracting <f>d from experimental data acquired under various conditions. Specifically, 
we are interested in the accuracy of these techniques over (i) the full range of differential 
phase cj) d E [0, 7 r], and (ii) a broad range of interferometer scale factor ratios k = Sk/SW 
To investigate these two aspects, we recorded data using the symmetric interferometer 
configuration shown in figure 4 with different interrogation times, T R b and Tk- Since k 
is proportional to (T K /T Rb ) 2 , each configuration of Tj corresponds to a different scale 
factor ratio. Additionally, the differential phase is modified with each Tk due to a 
systematic phase shift of the potassium interferometer from an external magnetic field. 
Therefore, we are able to study both effects with a single data set. 

Figure 8 shows a comparison between Bayesian and FRAC estimates of (fid, using 
the FRAC estimate as a reference. We varied T R b from 1 to 5 111 s, and Tk independently 
in the vicinity of T R b such that the scale factor ratio was modified over a relatively broad 
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range (k — 0.45 to 1.01). The phase noise due to the externally applied vibrations was 
kept quasi-common-mode between sensors by ensuring that Tk was within a few 100 /zs 
of Tr b- Over this range of Tr b and Tr, we found that the differential phase ranged from 
roughly <f>d = 0 to 2.8 rad as a result of a systematic shift of the potassium interferometer. 
It is clear from figure 8 that there is a high degree of correlation between the Bayesian 
and FRAC estimates, which is consistent with our expectations based on the simulations 
discussed in sec. 3.2. The error bars in this figure were computed from the combined 
statistical uncertainties of both methods, which both typically yield 8<f>d ~ 30 rnrad at 
the present level of noise. 

To summarize, we find that the difference between the two estimates is consistent 
with zero within a typical total uncertainty of ~ 40 rnrad. These data confirm that the 
two analysis techniques produce unbiased estimates of (f>d for dual-species interferometers 
with vastly different scale factors. We discuss further the advantages and limitations of 
these two techniques in the following section. 

6. Advantages and limitations of the methods 

As discussed in sec. 3.2, Bayes’ method is optimally efficient and yields a statistical 
error that scales as 1/\ f~N, compared to more heuristic fitting techniques which converge 
more slowly. This improved efficiency is a clear advantage of the Bayesian estimator 
compared to the FRAC analysis. However, the disadvantage is that it requires a priori 
information about the system, such as noise levels and interferometer contrasts, and 
it requires significant computational resources to evaluate. Furthermore, it is only a 
viable solution for simultaneous interferometer configurations that exhibit a high-degree 
of phase correlation. 

In contrast to the Bayesian estimator, the FRAC method requires only the 
interferometer timing parameters and a sensitive accelerometer that is well-coupled 
to the reference mirror in order to function accurately. It does not assume any 
particular interferometer configuration or require any additional system information. 
The FRAC method also has applications in absolute interferometry, as has been 
previously demonstrated in refs. [26, 66, 67, 68, 69]. Additionally, it is fast enough to be 
used for real-time feedback, which has been shown to improve single-sensor sensitivity 

[70]. 

Table 1 contains estimates of the phase noise for a single 8 'Rb interferometer, and 
two coupled interferometers of 8 'Rb and 39 K. The rms phase spread due to vibration 
noise (^) and the self-noise of the accelerometer (£) are shown for various frequency 
bands and interrogation times. For a single sensor analyzed with the FRAC method, the 
vibration-induced noise 4>fms represents the spread of phase on the uncorrected fringes, 
while the quantity (jx^f s indicates the residual phase noise present on the corrected 
fringes. Since this term is directly linked to the intrinsic noise of the mechanical 
accelerometer, it represents a fundamental limitation of the standard FRAC method. 
To give a quantitative example, based in the self-noise of the Titan accelerometer used 
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Single Sensor 



Differential Sensor 


Tkb (s) 

Noise 

Source 

DC - 

Frequency Range (Hz) 

1 1-10 10-100 DC - oo 

DC - 

Frequency Range (Hz) 

1 1-10 10-100 DC - oo 

0.01 

U vib ) 

0.023 

0.779 

1.071 

1.324 

0.000 

0.000 

0.004 

0.004 

U self ) 

0.000 

0.002 

0.002 

0.003 

0.000 

0.002 

0.002 

0.003 

0.1 

U vib ) 

Vrrms/ 

2.240 

19.62 

2.962 

19.97 

0.001 

0.380 

0.140 

0.405 

U self ) 

vrrms/ 

0.050 

0.078 

0.005 

0.093 

0.050 

0.078 

0.005 

0.093 


(cb vlh ) 
Vrrms/ 

58.13 

23.42 

3.779 

62.79 

1.067 

2.556 

1.511 

3.155 

1 

(<0 

2.936 

0.160 

0.006 

2.940 

2.936 

0.160 

0.006 

2.940 

(AT 

28.06 

1.868 

0.021 

28.12 

0.303 

0.098 

0.009 

0.319 


« f s )* 

0.294 

0.016 

0.001 

0.294 

0.294 

0.016 

0.001 

0.294 


Table 1. Comparison between phase noise for a single interferometer, and two coupled 
interferometers with effective wave vectors and . The rms phase noise (in 
radians) due to vibrations (i))*) and the self-noise of the mechanical accelerometer 
('(’rms) are shown for different frequency bands and interrogation times, Tr b- The noise 
from each band is summed in quadrature to obtain the total noise. Contributions less 
than 1 rnrad are not shown. For the simultaneous differential sensor, it is assumed 
that The rms phase noise <f%£ s was computed from eq. (C.9) using 

model (C.13) for the power spectral density S a (uj) of ground accelerations in a “quiet” 
location [66, 67] with integrated rms noise 1.4 x 10 -4 g [see figure C2(b)[. The quantity 
< t>s^ s was computed in a similar manner by replacing S a {w) with the self-noise spectrum 
of the accelerometer—here assumed to be white noise with |»S a | 1 / 2 ~ 3.2 x 10" 8 g/y/Hz. 
Quantities in the last row indicated by correspond to low-noise conditions that can 
be achieved with passive vibration isolation (integrated rms noise 1.4 x ICC 6 g ), and 
an accelerometer with 10 times smaller self-noise of ISal 1 / 2 — 3.2 x 1CT 9 g/y/Hz. 


in our experiments (3.2 x ICC 8 g/V Hz), the corresponding phase noise reaches ~ 90 
rnrad for an interrogation time of 100 ms, and ~ 3 rad by T kb = 1 s. With this level 
of self-noise, fringes cannot be reconstructed accurately. However, for a state-of-the-art 
device with an order of magnitude smaller self-noise (3.2 x 10~ 9 g/\J Hz), the phase noise 
decreases by a further factor of 10—allowing a least-squares fit to accurately converge on 
the fringe phase. We also point out that the noise contributions from both vibrations and 
self-noise are smallest at high frequencies—a result of the natural low-pass filtering of 
atom interferometers. Thus, high-bandwidth accelerometers are generally not required 
to implement the FRAC method with a single sensor. For instance, for ground-based 
WEP test facilities targeting a few ICC 15 [43, 44, 46] using interrogation times of T ~ Is, 
we estimate that a mechanical accelerometer with a self-noise less than ICC 11 g/y/Hz in 
the DC to 1 Hz frequency band will yield a phase noise contribution below the projected 
shot-noise limit of ~ 1 rnrad for each interferometer. 

When employing the differential FRAC method with two simultaneous 
interferometers, the self-noise of the accelerometer contributes to the phase of both 
sensors. Thus, in table 1 (where k = 1 and Tk = Tru yAfrb/fcjf ) > we have indicated 
the same values for in the corresponding columns for both single and differential 
sensors. However, we emphasize that one can measure the differential phase significantly 
more accurately than the self-noise limit for a single sensor. This is because the noise 
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introduced by the accelerometer is correlated between the two interferometers—reducing 
the uncertainty in the determination of cpd, as discussed in sec. 3.3. A recent study 
[71] has shown that uncertainties close to the quantum-projection-noise limit can be 
obtained with this method when the interferometers are in phase (<f>d = 0 ) and the 
accelerometer exhibits a conservative level of self-noise (0®ms ~ 0-3 rad). In general, 
for the differential FRAC method to function well for all values of <j>d and k, the self¬ 
noise of the accelerometer should correspond to less than n/2 in phase noise for each 
sensor—allowing individual fringes to be accurately fit. However, for the special case of 
4>d = 0 and k = 1 , the requirements on the accelerometer noise are much less stringent. 
For instance, ref. [71] indicates that reliable fits can be obtained with up to ~ 20 
rad. For these reasons, we emphasize that state-of-the-art mechanical accelerometers 
are not required to make sensitive measurements of (j)d with long-baseline differential 
interferometers. We anticipate that competitive levels of accuracy can be achieved with 
compact devices that feature a moderate level of sensitivity. 

For two coupled interferometers exhibiting different wave vectors, the vibration- 
induced phase noise is not identical and thus cannot be perfectly rejected at all 
frequencies. The values of listed in the last four columns of table 1 contribute 
directly to f>d —representing the level of uncorrelated differential phase noise in the 
system. We estimate that by = 1 s the differential phase noise reaches a level of 
~ 3 rad. However, we note that the differential transfer function [eq. (C.ll)] rejects 
most efficiently at frequencies below ~ 1/T, and this estimate is directly linked to 
the vibration spectrum used. In a quieter environment, such as that achieved with 
a vibration isolation platform [ 66 , 67] or in a satellite [52], the phase noise can be 
reduced by an order of magnitude or more. At this point, the Bayesian method can 
be employed—which easily handles differential phase noise. Since the sensitivity scales 
as (^rms/V^V, the analysis simply requires more measurements for larger to reach a 
given level of precision. 

7. Conclusion 

We have described and demonstrated experimentally two new analysis techniques for 
extracting the differential phase from coupled atom interferometers with different scale 
factors, Sj. A non-unity ratio n = S 1 /S 2 can result from using atoms with different 
kj, or from interferometers with different interrogation times, T y We also carried out 
correlated phase measurements between simultaneous interferometers of two chemical 
elements exhibiting different scale factors, and we have demonstrated a vibration 
rejection factor as large as 7 ~ 730. This system was used to validate the Bayesian and 
FRAC analysis methods, as well as a new ellipse fitting procedure [59], for extracting (pd- 
Furthermore, the FRAC method was used to demonstrate a statistical sensitivity for the 
Eotvos parameter of Sr) = 1.6 x 10 ~ 6 per measurement with T = 10 ms interferometers 
operating in a “noisy” environment. 

Presently, the contrast and SNR of our 39 K interferometer fringes are limited 
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by, respectively, the temperature of the atoms and technical noise present in the 
slave laser used for cooling, interferometry and detection. The latter issue results in 
uncorrelated offset noise, which reduces the correlation between interferometers. This, 
in turn, increases the error on cfd, which degrades the vibration rejection factor 7 
and the sensitivity Sr). We anticipate that a modest reduction of both the technical 
noise, and further cooling the 39 K sample in a gray molasses [79, 80], will result in a 
substantial improvement in the correlation between rubidium and potassium. A precise 
determination of 7 with our apparatus, including a complete evaluation of systematic 
effects, is beyond the scope of this work, but will be the subject of a future publication. 

Both the generalized Bayesian and differential FRAC methods yield unbiased 
estimates of cj)d for any scale factor ratio, k,, and are robust against experimental 
parameters such as the common phase range scanned by the two interferometers, or 
the level of uncorrelated offset noise present in the system. These features make both 
methods ideal for applications of dual-species interferometry where, until now, the 
available analysis tools could accommodate only systems that exhibit either k — 1 or low 
levels of common phase noise. These new methods are also appealing for gradiometer 
configurations using the same atoms and the same Tj [71], which have previously been 
utilized for precisely measuring the gravitational constant G and gravity gradients 
[27, 28, 29, 31, 32, 33, 61]. 

The freedom to vary the scale factor, the interrogation time or phase of either 
interferometer independently can be advantageous for studying systematic effects, 
interactions between atomic species [55], or for shifting the differential phase toward 
a region of higher sensitivity. Examples of such regions include <f>d — tt/ 2 in the 
case of ellipse-fitting methods, and f>d = 0 or n for the FRAC technique [71]. Both 
the FRAC and Bayesian methods also eliminate the systematic shift introduced on 
the measurement of A a when using dual-species interferometers with k 7 ^ 1—making 
them well-suited for upcoming WEP tests on ground [43, 44, 45, 46], in microgravity 
[41, 47, 48, 49, 50], and in Space [7, 51, 52, 53, 54], 
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Appendix A. Ellipse fitting methods 

In this appendix, we give some background regarding ellipse-fitting techniques and 
illustrate the problem of parameter bias for two different fitting algorithms. 

The general form of an ellipse in a cartesian plane is described by the algebraic 
equation for a conic 


F( A, y) = A • y = Ay\ + Bypy 2 + Cy\ + Vy 1 + £y 2 + J r = 0, (A.l) 


provided that B 2 < 4 AC. Here, A = {A, £>, C, V, £, J 7 } and y = {yf,yiy 2 ,y%,yi,y 2 A}- 
The center, orientation, major and minor axes of the ellipse are determined by the 
elements of A, and the differential phase can be shown to be 



(A.2) 


Generally, two types of ellipse-fitting algorithms exist: those that seek to minimize (i) 
an algebraic distance or (ii) a geometric / orthogonal distance between the ellipse and 
the data points. While algebraic methods tend to be simple, efficient and can guarantee 
an ellipse solution to the conic equation (A.l) [i.e. parabolic and hyperbolic solutions 
can be eliminated), they tend to suffer highly from bias in the ellipse parameters— 
resulting in a poor fit under certain circumstances. Geometric methods are usually 
much more accurate than algebraic algorithms, but at the cost of more complexity, 
more computation and less stability. Since minimizing the orthogonal distance between 
a point and an ellipse has no closed-form solution, these routines resort to iterative 
techniques that are not guaranteed to converge on an ellipse. 

A commonly used algebraic method is the simple and robust “direct ellipse fitting” 
(DEF) method developed by Fitzgibbon et al [62] that minimizes the sum of squared 
algebraic distances between the points and the ellipse, YliLi yf) 2 -, subject to the 
constraint B 2 = 4 AC — 1. Recently, Szpak et al [59, 63] developed an algorithm based 
on the optimization of the approximate maximum likelihood distance which seeks a 
balance between the costly geometric methods and stable algebraic techniques. This 
algorithm—termed the “fast guaranteed ellipse fitting” (FGEF) method—also includes 
error estimation for the geometrically meaningful ellipse parameters (center coordinates, 
axes and orientation) which we have extended to include an estimate of the differential 
phase error, 8<f>d- 

Figure Al illustrates the bias introduced on the differential phase estimated by the 
DEF and FGEF methods. For moderate amounts of noise in the offset, the DEF method 
tends to produce fits that are characteristically compressed along the major axis and 
stretched along the minor axis of the ellipse, as shown by the red curve in figure Al(a). 
This effect results in a biased estimate of (fid that increases monotonically away from 
7 t/2 , as shown in figure Al(b). In contrast to the DEF method, the FGEF algorithm 
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Figure Al. (a) Synthetic data following an ellipse with added offset noise. The solid 
green curve represents the actual ellipse, and fits to the data using the DEF method 
(red curve with big dashes) and FGEF method (blue curve with small dashes). The 
simulated ellipse contains 500 points with Gaussian-distributed noise on the offset 
parameters Bj with standard deviations {<JB 1 ,<rB 2 } = {0.01,0.03} (corresponding to 
SNR ~ {20,6}). Ellipse parameters: A x = A 2 = 0.2, Bi = B 2 = 0.5, k = 1, 
fd= 1 rad. (b) Measured bias in differential phase estimates, (jff 1 , from the DEF (red 
triangles) and FGEF (blue points) methods relative to the actual value, (j> d ct . The 
black squares show the estimates from the differential FRAC method for comparison. 
On all plots, the error bars correspond to the statistical distribution of fits to 100 
synthetic data sets. 


predicts an ellipse (shown in blue) that is much more representative of the actual ellipse 
(shown in green), and also results in less bias in (j) d in the central region around 7t/2. 
Outside of this region, the bias behaves nonlinearly in a manner that depends on the 
ellipse parameters and the level of noise. Here, we point out that these bias estimates 
are dependent on the type of noise (offset, amplitude, or differential phase) and the 
amount of noise present in the data, but typically the bias is smallest in the vicinity 
of (j) d = 7 t / 2 , and decreases with the noise level. In general, ellipse-fitting techniques 
always generate a non-zero systematic on the differential phase estimate, and depending 
on the level of sensitivity, this bias must be carefully accounted for when performing 
precise measurements with f> d [28, 31, 32, 34, 61]. 

Appendix B. Bayesian analysis of Lissajous curves 

In this appendix, we describe in detail our generalized Bayesian analysis technique 
to estimate the differential phase from Lissajous curves. We also demonstrate the 
effectiveness of this method using numerically simulated data with Gaussian noise in 
the offset parameters {Bi,B 2 } and the differential phase, <f d . Noise in the amplitude 
parameters {Ai : A 2 } of the coupled-sensor model (2) can also be included via a trivial 
modification of the noise model. In what follows, we first provide some relevant 
theoretical background of the Bayesian estimation technique. For a more comprehensive 
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description of Bayesian analysis in this context, see ref. [65]. 

In a generalized system, where M represents a measurement of the system quantities 
and V represents a variable we are interested in measuring, Bayes’ rule can be 
summarized by the following equation 


P(V\M) 


p(V)L(M\V) 

N(M ) 


(B.l) 


Here, P(V\M) is called the “posterior” probability distribution and represents our 
state of knowledge after a measurement, M. p(V) is the “prior” probability before 
the measurement, and L(M\V) is called the “likelihood” to obtain a certain result for 
M given V. The key to the entire estimation process is the likelihood distribution, 
which is computed based on a specific model of the noise present in the system. The 
quantity N(M) = ^2 V L(M\V)p(V) is the probability of measuring M integrated 
over all possible values of V, and is just a normalizing factor for the posterior 
distribution. Mathematically, L(M\V) can be thought of as a function of V with M 
fixed, and vice versa for P(V\M). The essence of Bayes’ rule is that knowledge of 
the variable V can be updated on a measurement-by-measurement basis—with each 
successive measurement contributing additional information that narrows the width of 
the probability distribution associated with V. A well-known example of this type 
of recursive analysis is a Kalman filter [81], which is used extensively in the fields of 
guidance, navigation and trajectory optimization. 

For the specific case of two coupled atom interferometers, the variable of interest 
is <f>d and the i th system measurement is given by the pair of (normalized) atomic state 
populations Mi = {ni,n 2 }i- Thus, for a single measurement eq. (B.l) becomes 


P{M{ n ^n 2 } i ) i = 


p((j)d)iL({n u n 2 }i\(j)d) i 


(B.2) 


N({n 1 ,n 2 }i) i 

where P(f>d\{ni,n 2 }i)i is referred to as the conditional distribution based on the i th 
measurement. The basic algorithm for Bayes’ estimation can be summarized as follows: 


(1) Choose a suitable initial prior distribution, p(<f>d)i= i- In our case, we take this to 
be a uniform distribution within the range f>d £ [0,7r], and zero elsewhere. 

(2) Record a new measurement {ni,n 2 }i, and calculate the likelihood distribution 
L({rii, n 2 }i\(pd)i from the noise model. 

(3) Compute the conditional probability distribution P((f>d,\{ni, n 2 }i)i from Bayes’ rule 
(B.2). 

(4) Set the new prior distribution equal to the previous conditional distribution: 
p(4>d)i +i = T , (0 d |{?zi,n 2 }i) i . 

(5) Repeat steps (2) through (4) until the width of the conditional distribution reduces 
to the desired level. 


(6) Estimate the variable of interest, (fd, using the maximum likelihood value of the 
final conditional [i.e. posterior) probability distribution. 

This algorithm is illustrated in figure 1 in sec. 3.2. 
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The main challenge in Bayesian analysis is to compute the likelihood distribution 
L({ni, 712 }|0d) given a specific model for n± and 77 . 2 . For the specific case of coupled 
interferometers, there are three possible sources of noise: amplitude, offset and 
differential phase. To illustrate each source, we modify the definitions of the nj in 
eq. (10) to explicitly include these noise terms 

ni(4> c ) = (1 + Mi) cos (k4> c + 4> d + tyd) + fiBi, (B.3a) 

n 2 (f)c) = (1 + SA 2 ) cos( 0 c ) + SB 2 . (B.3b) 


The parameters fiAj, fiBj, and fif>d represent uncorrelated noise in the amplitude, offset 
and differential phase, respectively, each of which is assumed to follow a Gaussian 
probability distribution with zero mean and non-zero standard deviation. Using this 
model, the likelihood distribution can be shown to be [65] 


L({7i 1 ,n 2 }\(pd) 


f 1 P(ni\si)P(n2\{s2/,Si,f>d}) 

i J-i V 1 - s i 


(B.4) 


Here, P(n 1 |s 1 ) and P{n 2 \{s 2 /, Si, f>d}) are the single-sensor conditional probability 
distributions for ri\ and 77,2, which we discuss in more detail below. The quantities 
Si = cos(k0 c + <f>d ) and S2,i = cos(0 c ) are the principle variables on which the coupled 
measurements 77,1 and 77,2 depend in the model (B.3). Due to the periodic nature of 
the Lissajous equations (10), for each value of ri\ there are multiple possible solutions 
for n 2 (as shown in figure Bl). We assign an integer i to each of these solutions. 
More specifically, S 2 / is the f th root of 112 given n\ = sq. The sum over l appearing 
in eq. (B.4) accounts for all possible solutions. In the distribution functions P(ni|si) 
and P(n 2 \{s 2 /, Si, f>d}), we denote the implicit dependence on variables Si and f> d by a 
semi-colon. This notation emphasizes that the quantity S 2 / is coupled to Si through the 
common phase f> c ^[. Finally, we point out that the coupled variables Si and s 2 / both 
depend on <pd , but we do not write this dependence explicitly. 

At this point, we need to know the possible values 7i 2 = S 2 / (given a measurement 
of n-i = Si) which enter into the likelihood distribution. We devote the remainder of this 
section to a detailed description of computing the roots of the Lissajous equations (10). 
As mentioned above, due to the non-linear nature of Lissajous curves, there are multiple 
possible solutions for 712 given a single value of n\ within a predefined phase range. We 
denote these solutions s 2 ,£ for integer L When k — 1, the Lissajous curve collapses to 
an ellipse, and only two values of 772 exist for each 77i over any 2n range of f> c . In this 
case, it is straightforward to compute the two solutions as s 2 ,±i = cos[cos^ 1 (si) ± <f>d\. 
However, when k 7 ^ 1, the problem is much more complex. If the scale factor ratio can 
be written in the form k = p/q, where p and q are prime numbers, then the period of 

Since it is assumed that 4> c is random and unknown, the probability distributions of Si and S 2 / are 
equivalent to that of a sinusoid: -P(s2,dbc) = (1 — s\ ■ 




Correlative methods for quantum tests of the weak equivalence principle 


29 




n\ 


n i 


Figure Bl. Example of a Lissajous curve for k = p/q = 3/7. (a) The curve is plotted 
over the phase range 4> c £ [0, 27t] and the solutions { ni , n 2 = s 2 ,e} are shown as points. 
For n\ = —0.9 there are two solutions s 2 /, shown as the blue points. Similarly, for 
ni = 0 there is only one possible value of n 2 (shown in green), and for n\ = 0.8 no 
solutions exist, (b) The same Lissajous curve plotted over (j) c £ [0, 27 t< 7 ] with q = 7. In 
this range, there are always 6 solutions S 2 ,i for each value of n\ (although they may 
not be unique). 


the Lissajous curve is 2nq —requiring q revolutions to form a closed loop. Within each 
27r interval, there can be either 0, 1 or 2 solutions of n -2 for each n \, as illustrated in 
figure Bl. 

To calculate these solutions for a given ri\ = si and 4>d, it is necessary to know the 
approximate range of common phase spanned by the data: 0 C G [</>™ m , 0™ ax ] + . With this 
information, we compute the range of phase spanned by sensor 1, 9 G /t[</>“ in , 0 aiax ] + 4>d, 
and we subdivide this range into intervals of n such that the £ th interval is defined as 
the range G [l,d + l)7r, where £ = [6 1 /7rJ. Here, the brackets |_- • ■ J indicate the 
floor function. Beginning with the left-most interval, we check for solutions sequentially 
at each n phase bin until the entire range is spanned. Empirically, we find that if a 
solution exists within the K th interval given a value n\, then it is unique and can be 
written explicitly as 


cos[(cos 1 si — (fd + 2nmi t e)/K\ for even £, 

' cos[(cos _1 si + <pd — 27T m 2 /)/n\ for odd £, 

where the integers m\^ and rri 2 j are defined as 

L(£ - i)/2j e<-i, ( K/2J e<-2, 

m\£ — \ 0 — 1 < £ < 2 , m 2) i —l 0 — 2 < £ < 1 , 

[£/2\ £>2, [(£ + l)/2\ £>l. 


(B.5) 


(B.6) 


+ In practice, this range estimate does need not to be very precise—we find that estimating the correct 
range to within ±7r still results in a precise estimate for fd- However, overestimating the phase range 
may result in a slower convergence rate for the estimate. See sec. 5.1 for a description of how the phase 
range can be estimated experimentally. 
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Figure B2. Simulated Bayesian phase estimates in the presence of offset noise. Top 
row: simulated data for re = {0.6, 0.8,1.0,1.2,1.4}. Here, points based on eq. (B.3) 
were generated from common phases chosen randomly over the range 4> c € [—27t, 27t] . 
The differential phase was set to (f>^ ct = 1 rad, and offset noise was applied to the 
points via Gaussian distributions with standard deviations {c tb 1 , ctb 2 } = {0.02,0.04} 
corresponding to SNR ~ 50 and 25, respectively. Bottom row: the systematic error in 
the Bayesian phase estimate, (f >\} st — as a function of for each value of re. Only 
10 points were required to bring the statistical uncertainty indicated by the error bars 
to < 50 mrad. 


With these solutions in hand, it is possible to compute the likelihood (B.4) given 
specific noise models for the single-sensor probability distributions P(ni|si) and 
P(n 2 \{s 2 /] s i, cpd})- We now investigate the specific cases of offset and differential phase 
noise on the extraction of fd from simulated data sets. This analysis can also be extended 
to include noise in the fringe amplitudes through the parameters 5Aj [65], but we do 
not consider this case here. 

Offset noise 

When the system exhibits noise only in the offset of the atomic state measurements, 
the parameters SBj are randomly distributed for each repetition of the experiment, 
and SAj = 0 and 5<f>d = 0 in the model (B.3). Under realistic conditions, these noise 
parameters follow a Gaussian distribution with zero mean and standard deviations given 
by erg., and the single-sensor conditional probabilities can be written as 

P(ni\ Sl ) oc exp [ - (m - s l ) 2 /2a 2 Bi ] , (B.7a) 

P[n 2 \{s 2 fi si, fd}) oc exp [ - (n 2 - s 2 y) 2 / 2 ^I 2 j • (B.7b) 

Figure B2 shows some examples of simulated data in the presence of offset noise, 
where the differential phase has been extracted using the Bayesian estimation algorithm 
described above. These simulations show that <f>d can be precisely estimated over the 
full range of 0 — 7r, and for a wide variety of scale factor ratios. Here, we demonstrate 
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Figure B3. Simulated Bayesian phase estimates in the presence of differential phase 
noise. Top row: simulated data for k = {0.6,0.8,1.0,1.2,1.4}. Similar to figure B2, 
common phases are randomly chosen over the range <j> c £ [ — 7r, 7t] , and the differential 
phase is set to <p d ct = 1 rad. Differential phase noise is applied to each point with a 
standard deviation of a ( j >d =0.1 rad (SNR ~ 10). Bottom row: the systematic error 
in the Bayesian phase estimate, q !>^ st — 4> d ct , as a function of (f> d ct for each value of k. 
Only 10 points were required to bring the statistical uncertainty to < 50 mrad. 


the technique for the limited range k G [0.6,1.4], but we have also verified that the 
extraction method works well outside this range. In contrast to ellipse-fitting techniques, 
no systematic bias in the phase estimates is observed, and fewer points are required to 
converge to competitive error levels. 

Differential phase noise 

Since the noise parameter associated with the differential phase, Scpd, adds directly to 
the quantity of interest, fa, we can account for this type of noise by adding an extra 
convolution with our noise model at the end of any likelihood calculation. We choose to 
examine the case of Gaussian noise for the differential phase, such that the conditional 
probability distribution is 


P^'dl^d) oc exp[-(0' d - (f d ) 2 /2al d }. (B.8) 

Here, (f)' d represents a measured value of the differential phase in the presence of Gaussian 
noise centered on the most likely value, f>d, and a r p d is the standard deviation of the noise 
distribution. The modified likelihood function is described by the convolution 

/ OO 

P({ n l: n 2}\(p'd) Pi^dl^df^'d- (B.9) 

-OO 

In a similar fashion to the offset, in the absence of any other noise sources it is necessary 
to estimate multiple candidate solutions for (fd over a given range of <f c in order to 
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compute the likelihood function. Before convolving with the conditional probability 
distribution in eq. (B.9), the likelihood function can be written as 

L({n u n 2 }\f' d ) = - <f>d,k), (B.10) 

k 

where S(x) is the Dirac delta function, and the sum over k accounts for all candidate 
solutions (j>d,k that exist in the common phase range <p c G [0™ m , 0™ ax ]. These solutions 
can be computed by, again, dividing the phase range into intervals of 7r, and labeling 
each of them by an integer k = |_0 C /7rJ. We find that two possible solutions exist for 
<; t>d within each interval, which we denote as (p'^f for (p c G [±kn, ±(fc + l)7r). Explicitly, 
these phases can be computed from 

= cos _1 (rii) ± K[cos _1 (n 2 ) - 2nm k \, (B.ll) 

where = (—l) fc |_(| k\ + 1)/2J. We transform these phases into the range of 0 — 7r using 
^tk cos _1 [cos(0^)]. Finally, this result is convolved with the Gaussian noise model 
to obtain 

L({ni,n 2 }\(j) d ) = ^exp[-(0 d - ^J) 2 / 2 ^]- ( B - 12 ) 

k 

Two subtleties exist with this analysis, however, that warrant discussion. First, 
when the common phase range exceeds (p c G [—7r,7r], the Bayesian analysis may predict 
multiple equally probable values for (pd- This is obviously a problem if we are interested 
in a precise, unique estimate of the differential phase, and we have no pre-existing 
knowledge of its value. Therefore, we restrict our consideration of the problem to a range 
of common phase within — 7 r to 7 r. Second, the noise parameter 5(pd can theoretically take 
any value, i.e. 8<pd G (— 00 , 00 ), although in practice it is limited to a finite range defined 
by <J r j, d . So far, we have considered (pd only in the range of 0 to 7 r, but for situations 
where a ( j >d > 7 t/4 , the likelihood distribution can have significant contributions from 
the wings of the adjacent 7r phase intervals. This effect can be taken into account by 
using the fact that P(<pd) = P{—(pd) = P(2 tt — (pd), and adding mirrored versions of the 
likelihood to the convolution in eq. (B.9). This “tiling” technique can be extended to 
account for large noise levels, where more than one 7 r phase bin is spanned [65]. 

Figure B3 shows some examples of simulated data in the presence of differential 
phase noise. As for the case of offset noise, estimates of <pd exhibit no significant bias 
over the full range of 0 — 7r, and for a large range of scale factor ratios. Additionally, 
only a small number of points are required to converge to a level of uncertainty less than 
that of the noise defined by a ( f >d . The convergence of this uncertainty as a function of 
the number of measurements is the subject of the next section. 

Scaling with measurement number 

To test the scaling of the statistical and systematic error of the Bayesian estimator 
as a function of the number of measurements, we performed the following study. We 
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Figure B4. The statistical and systematic error in the Bayesian phase estimator as 
a function of the number of measurements, N, used in the Bayesian analysis for (a) 
differential phase noise and (b) offset noise. In both plots, k = 0.8, <p^ ct = 1 rad, 
and M = 50 samples were used. Black points represent the statistical uncertainty 
, and blue triangles indicate the systematic error ejf. The solid red lines indicate 
the minimum convergence rates based on a^ d l\jN for differential phase noise and 
on 1 / \JNIIfd) for offset noise, where I(4>d) is the Fisher information given by 
eq. (B.14). The dashed horizontal lines represent the nominal phase resolution used in 
the simulations. Values of cr c / )d = 0.1 rad (SNR ~ 10) and cjb-l = Tb 2 = 0.2 (SNR ~ 5) 
were used as the noise parameters in (a) and (b), respectively. 




randomly generated M = 50 samples of “measurements”, each containing 100 points 
following the model (B.3) with noise added to either the differential phase or the offset. 
As a function of the measurement number, N, within each sample, we computed the 
Bayesian estimate 0^ st (Af) and the standard deviation of the associated probability 
distribution The statistical error for each measurement is taken as the average 

of S(f^ t (N) over all M samples, which we denote as e^ d at (lV) = (Scf> e c f t (N)) m- Similarly, 
the systematic error is defined as e s J^(N) = (|</>)) st (lV) — <^ ct | )m- The results are shown 
in figure B4. 

For the specific case of noise that contributes directly to the variable of interest 
{e.g. differential phase noise) the statistical uncertainty of the Bayesian estimator is given 
by e^ at = (jfa/y/N. As we show in figure B4(a), the measured statistical error closely 
follows this dependence. Similarly, on average the systematic error drops to a level much 
less than e^ tat after only a few measurements. This level is primarily determined by 
the grid resolution used when computing the likelihood distribution for (j)^. During the 
estimation procedure, we initially set the phase grid resolution to ~ 7r/ 100, and we refine 
this grid size on a measurement-by-measurement basis. As the likelihood distribution 
narrows, grid points are redistributed toward the maximum likelihood value. We find 
that this grid optimization procedure can improve the resolution by up to an order of 
magnitude (depending on the level of noise in the system), while keeping the number of 
integral evaluations per measurement fixed. 

For the more general case of noise present in a parameter that is indirectly related 
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to the quantity of interest through some function, the uncertainty is constrained by the 
Cramer-Rao lower bound ^ 

(B-13) 


stat \ 

e <t> d - 


y/NIfa) 


This relationship can be used to compute the minimum convergence rate of the statistical 
error in the presence of offset or amplitude noise, for example, where the noise affects 
4>d indirectly through the quantities { 711 , 712 }. The Cramer-Rao lower bound includes 
the Fisher information, I(<f>d), of an individual measurement, which can be computed 
from the likelihood distribution L({ni,n 2 }\(fd) as follows 




d 2 

W d 


111 [L({?7i,77 2 }|0d)] 


{ni,n 2 } 


f f' dn 1 dn 2 

( dL ) 

JJ L({m,n 2 }\(f> d ) 

\d(pd) 


(B.14) 


Here, the brackets (• • • ){n 1 ,n 2 \ denote an average over the random variables {n 1 ,7i 2 }- 
The Fisher information is a measure of the amount of information that a random 
variable (or a set of random variables) carries about an unknown parameter. In this 
case, the unknown parameter of interest is (j) d and the set of random variables is 
the set of measurements {tii, 71 2 }, which are governed by the likelihood distribution 
L({ni, n 2 }\(j)d) —hence its appearance in eq. (B.14). This quantity has no closed-form 
expression for the case of offset or amplitude noise in our system, and must be evaluated 
numerically. For the parameters used in figure B4(b), we find I((j> d ) — 3.3, which gives a 
minimum convergence rate of 0.55 /\/N. This rate is consistent with the measured 
statistical uncertainties shown in the figure. We note that the Fisher information 
empirically scales as I ~ e ~ /3aB , where ft is a large factor that depends on the differential 
phase and the scale factor ratio used (e.g. ft ~ 35 for n = 0.8 and <f> d = 1 rad). Thus, with 
only a moderate reduction to the level of offset noise in the system, one can dramatically 
improve the convergence rate of the Bayesian estimate. 


Appendix C. Response of a dual-species interferometer to mirror vibrations 


Here, we summarize the essential theoretical tools required to evaluate the response of 
both single- and dual-species interferometers to vibrational noise of the retro-reflection 
mirror. 

First, we provide a review of the sensitivity function for a single atom interferometer, 
g(t). This function characterizes how the interferometer transition probability behaves 
in the presence of fluctuations in the Raman laser phase difference, ip-hit). Developed 
previously for use with atomic clocks [82], the sensitivity function is a useful tool that 
can be applied, for example, to evaluate the response of the interferometer to laser phase 
noise [57], or to correct for spurious vibrations in the Raman beam optics [41, 58, 69]. 
We are primarily interested in the latter. 

The sensitivity function is a unitless quantity that is defined as follows 


9(t) 


lim 

8 ip —>-0 


5$(6ip, t ) 
dip 


2 lim 

>-0 


8 P(8ip, t ) 


Sip 


(C.l) 
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where Sip is a phase jump occurring at time t during the interferometer that modifies 
the total interferometer phase, <f>, by an amount <5$, and the transition probability 
P(<f>) = (1 — cos<h)/2 by a corresponding amount 5P. Thus, the interferometer phase 
due to an arbitrary phase noise function, <p(t), can be computed as 

(C.2) 

The quantum mechanical nature of the atom plays a crucial role on the sensitivity 
function—in particular, the evolution of the internal atomic states during each Raman 
pulse. Using the procedure outlined in refs. [41, 57], the sensitivity function, gj(t), of an 
interferometer with timing parameters labeled with subscript “j” can be shown to be 


9j{t) 


-sm(Qf (t- AT})) 

-1 

— sin — ATj — 7})) 

1 

- sin (Qf(t - ATj - 27})) 


0 < t — ATj < Tj, 

Tj < t — ATj < Tj + Tj, 

Tj + T j < l ~ ATj < Tj + 3rj, 
Tj + 3 Tj < t — ATj < 27} + 3tj, 
2 Tj + 3tj < t — ATj < 2 Tj + 4 Tj, 
otherwise. 


(C.3) 


ffere, Tj is the interrogation time, Tj is a pulse duration, is the effective Rabi 
frequency associated with the two-photon Raman transitions, and ATj is a delay with 
respect to t — 0 that facilitates a difference in the start time between interferometers. 
It is assumed that = n/2, such that the first and third interferometer pulses have 

pulse areas of 7 t/ 2 with duration Tj, and the second is a 7r-pulse of duration 2 Tj. 

To evaluate the response of an interferometer to Raman mirror motion, the phase 
noise function is first expressed as <Pj{t) = kjz(t ), with z(t ) representing the time- 
dependent position of the mirror along the axis of the beams. Then, the phase shift of 
interferometer j due to movement of the Raman mirror is 

/ OO poo 

Wj(t)a vih (t)dt = kf / fj(t)a vih (t)dt, (C.4) 

-OO J — OO 


where a vlb (t) = z(t ) is the time-dependent acceleration of the mirror due to vibrations, 
Wj{t) = kj S fj(t) is a time-dependent weight function for the mirror accelerations, and 
fj(t) is called the response function associated with the j th interferometer. This function 
is given by the integral of the sensitivity function: fj(t ) = — g ] {t')dt', and can be 
evaluated as 


fj (t) 


Aft (1 - cos Qf(t- ATj)) 

(t — ATj) + Aft — Tj 
Tj + j^(l - cosfif (t - A Tj - T)) 
27} + 3 Tj - — {I — ATj) 

^(l-cosnf(f-A7}-27})) 

0 


0 < t — ATj < Tj, 

Tj < t — ATj <T + Tj, 

Tj + Tj < t — ATj < Tj + 3 Tj, 
Tj + 3 Tj < t — ATj < 27} + 3 Tj, 
27} + 3tj < t — ATj < 2Tj + 4 Tj, 
otherwise. 


(C.5) 
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Figure Cl. Weight functions, Wj(t), described by the response function (C.5). These 
weights determine the phase shift associated with mirror vibrations in eq. (C.4). The 
pulse durations, Tj , satisfy = 7t/2. The differential weight function, i.e. the 

difference between the red and blue curves, is shown in black. 


At its heart, eq. (C.4) is a generalization of the well-known interferometer phase shift 
due to a constant acceleration, a 

0. = S ja = kf{Tj + 2 T j) (t, + a — kfTf a. (C.6) 

In this relation, the quantity Sj ~ is equivalent to the integral of the weight 

function, Wj(t), which determines how strongly the mirror vibration at time t contributes 
to the interferometer phase shift. This function is triangle-shaped, as shown in 
figure Cl, which indicates that the phase contributions are smallest near t = AT) and 
A Tj + 2 Tj + 4tj , where the wavepacket separation is a minimum. Similarly, the weights 
are largest near the mid-point, t = AT) + Tj + 2 Tj, where the separation between the 
interfering states is a maximum. 

For the case of two coupled interferometers, the differential phase shift resulting 
from mirror vibrations can be expressed as 

/ OO 

w d (t K»>(()di, (C.7) 

-oo 

where the differential weight function, w d (t), is given by the difference between the 
single-sensor weight functions 

w d (t) = wi(t) - w 2 (t ) = kffxit) - kff 2 (t). (C.8) 

This function has an intuitive understanding. For the extreme case when = /c| ff and 
the two interferometers are perfectly overlapped (i.e. A T x = AT 2 , T x = T 2 , r x = r 2 ), 
Wd(t) is zero everywhere. This implies that the differential phase shift due to mirror 
motion is 0^ lb = 0 —corresponding to perfect common-mode phase noise rejection. In the 
opposite extreme, when either ^ klf or the interferometers are not well-overlapped, 
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vibration noise induces a differential phase shift </>^ lb between the two sensors given by 
eq. (C.7). This non-zero phase shift is directly responsible for uncorrelated contributions 
to f>d in the case of non-overlapped interferometers, and it explains the loss of common¬ 
mode rejection in the case of coupled interferometer with different scale factors. For the 
case of a constant acceleration, eq. (C.7) can also be used to derive the systematic shift 
Sf> s J b = (Si — S 2 )a resulting from interferometers exhibiting S± ^ S 2 . 

One can characterize how mirror vibrations with a given frequency spectrum affect 
each interferometer by computing the mean-squared phase noise 

1 7°° 

(Cn b s) 2 = ^ J Q \Hj(uj)\ 2 S a (uj)duj. (C.9) 

Here, S a (uu) is the power spectral density of acceleration noise on the mirror, and Hj(uj ) 
is the transfer function associated with interferometer j given by the Fourier transform 
of Wj(t ). The transfer function describes how acceleration noise at a given frequency 
affects the phase over the duration of the interferometer. For frequencies c a <C and 
pulse separations Tj Tj, this function is well-approximated by 

Hj(u) = — ie-to&Tj+Tj+^kf Jf sine 2 (C.10) 

For the dual-species interferometer, one uses the differential transfer function in the 
same fashion: 

H d (u) = Hi(u) ~ H 2 (u). (C.ll) 

These functions are shown in figure C2(a) for realistic experimental parameters 
associated with a K-Rb interferometer. Here, there is a clear difference between the 
transfer functions associated with single-species and dual-species interferometers. For 
the individual sensors, the transfer function is well-approximated by the square of a 
sine function—which exhibits regular zeroes at the fundamental frequency 1 /Tj and an 
envelope that decreases as (2 /coT )) 2 . This dependence implies that the interferometer 
naturally filters the high-frequency components of the vibration spectrum, with a —3 
dB cut-off frequency of (n™ t /27r = \/2 /ttTj ~ 1/2 Tj. 

The differential transfer function, on the other hand, has a much more complicated 
frequency dependence. We will focus on the most interesting case for WEP tests, 
i.e. when the wave vectors satisfy = (1 — e) 2 fcf ff where e -C 1, and the two 
interferometers are symmetrically overlapped in time as shown in figure Cl. Under 
these conditions, we find that Hd can be approximated by 

\H d (uj)\ » kf(T 2 - T 2 2 )sinc ^ 

- 2ekf T 2 sinc 2 

It follows that there is a competition between the two terms in this expression. For 
the extreme case when e = 0 (i.e. kf s = ktf), the differential transfer function is 


oj(Ti + T 2 ) 


sine 


w(7i - T 2 ) 


(C.12) 
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(b) 



Figure C2. (a) Normalized transfer functions, \Hj(uj)\/kj S Tf, described by eq. (C.10) 
for coupled K-Rb interferometers. These functions determine the response of the 
single-species (red curve) and dual-species (black curve) interferometers to acceleration 
noise at different frequencies, to. Here, T) = 1 s and T 2 = (1 — e)Ti, with 
e = 1 — y^kfjkf — 0.0087, and n = t 2 = 10 /rs. (b) Model curves for the 
power spectral density of ground vibrations described by eq. (C.13). The solid red 
curve corresponds to the power spectral density of a “quiet” location [66, 67] that we 
simulated with the parameters £ = 0.5, % = 1, a>o/2i r = 10 Hz, a rms = 1.4 x 10~ 4 g. 
The dashed blue curve corresponds to “low-noise” conditions achievable with passive 
vibration isolation (£ = 2, x = 0-8, u}q/2tt = 1 Hz, a rms = 1.4 x 10 -6 g). These two 
curves were used to compute the rms phase noise in table 1. 


dominated by the first term, which is identically zero for all frequencies only if Tf = T 2 . 
This represents the ideal case for gravity gradiometry applications. On the other 
hand, when e > 0 and 7\ = T 2 , the second term in eq. (C.12) dominates. Since 
the two interferometers are assumed to have different wave vectors, it is not possible 
to make the transfer function zero at all frequencies. However, it is straightforward 
to show that Hd = 0 in DC provided that k^Tf = kf^Tf. This criteria optimizes 
the rejection of common-mode vibration noise at frequencies below the cut-off for a 
single-sensor, uij ut , and can be achieved by adjusting the interrogation times such that 
To = ^kfJWTi = (1 - e)Ti. 

Figure C2(a) shows a comparison between single-sensor and differential transfer 
functions for T ~ Is interferometers. When operated differentially, the sensitivity to 
vibrations at frequencies less than ojj Ut is typically more than 3 orders of magnitude 
below that of the single interferometer, despite the fact that kf s ^ k?r ■ 

Figure C2(b) displays the power spectral density function, S a (u>), that was used to 
compute the rms phase noise in table 1. These curves are based on a regression model 
for ground accelerations [83] described by 

S a (u>) = N(£,x) ^ 

This model has a single peak at cu = wo and contains two positive shape parameters, 
£ and y, that determine the sharpness of the peak and the scaling of the wings of 
the distribution, respectively. The quantity a^ ms is the mean-squared acceleration of 




u; 2 - UJ%) 2 + (2£o;ot<;) 2 
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the corresponding time-domain acceleration signal, aft), and N(£,x) is a normalization 
factor that depends on the shape parameters. This factor is chosen such that the integral 
f-oc = 2tt a r 2 ms . 
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